      Subroutine EULER
      include "param.fi"
      include "common.fi"

      Integer I

      !If (Shepard .EQ. 1) Then
        Call Calc_ShepardCoef
      !Endif

      If (ifDual.EQ.1) Then
        Do I = 1, GridDual2
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       Rho, RhoDual(I))
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U1, U1Dual(I))
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       U3, U3Dual(I))
        Enddo
        Do I = 1, NTOT3T
          Call IntepGrad(X1(I), X3(I), U1Dual, DU1DX(I), DU1DZ(I))
          Call IntepGrad(X1(I), X3(I), U3Dual, DU3DX(I), DU3DZ(I))
        Enddo
      Endif

      If (ifDual.EQ.1) Then
        Call Calc_dRho_dT_Dual
      Else
        Call Calc_dRho_dT
      Endif
      Do I = 1, NTOT3T
        Rho(I) = Rho(I) + Temp1(I) * DT
        If (Rho(I).LT.RDENSITYMIN) Rho(I) = RDENSITYMIN
        If (Rho(I).GT.RDENSITYMAX) Rho(I) = RDENSITYMAX
      Enddo

      !Call AVER(Rho, Temp1)
      !Do I = 1, NTOT3T
      !  Rho(I) = Temp1(I)
      !Enddo

      !Write(*,*) "Calc_dRho_dT", U1(1)

      Call PressureField
      !Write(*,*) "PressureField", U1(1)
      If (ifDual.EQ.1) Then
        Do I = 1, GridDual2
          Call Interpolation(X1Dual(I), 0.0, X3Dual(I),
     &                       P, PDual(I))
        Enddo
        Do I = 1, NTOT3T
          Call IntepGrad(X1(I), X3(I), PDual, DPDX(I), DPDZ(I))
        Enddo
      Endif

      Do I = 1, NTOT3T !- NPWALL
        ! Gravity
        U3(I) = U3(I) + Gravity * DT
      Enddo
      If (ifDual.EQ.1) Then
        Call Calc_dU_dT_Dual
      Else
        Call Calc_dU_dT
      Endif
      Do I = 1, NTOT3T !- NPWall
        U1(I) = U1(I) + Temp1(I) * DT
        If (DIMM.EQ.DIMM3D) U2(I) = U2(I) + Temp2(I) * DT
        U3(I) = U3(I) + Temp3(I) * DT
      Enddo
      !Write(*,*) "Calc_dU_dT", U1(1)

      If (ALE .EQ. ALEEuler) Then
        Call Grad(U1, DU1DX, DU1DY, DU1DZ)
        If (DIMM.EQ.DIMM3D) Call Grad(U2, DU2DX, DU2DY, DU2DZ)
        Call Grad(U3, DU3DX, DU3DY, DU3DZ)
        Call Calc_NLTerm
        Do I = 1, NTOT3T - NPWall
          U1(I) = U1(I) + Temp1(I) * DT
          If (DIMM.EQ.DIMM3D) U2(I) = U2(I) + Temp2(I) * DT
          U3(I) = U3(I) + Temp3(I) * DT
        Enddo
      Endif

      If (NU .NE. 0.0) Then
      Call Calc_VisTerm
      Do I = 1, NTOT3T !- NPWall
        U1(I) = U1(I) + Temp1(I) * DT
        If (DIMM.EQ.DIMM3D) U2(I) = U2(I) + Temp2(I) * DT
        U3(I) = U3(I) + Temp3(I) * DT
      Enddo
      !Write(*,*) "VisTerm", U1(1)
      Endif

      If (ALE .EQ. ALELagrange) Then
      Do I = NTOT3T - NPWALL + 1, NTOT3T
        U1(I) = 0.0
        U2(I) = 0.0
        U3(I) = 0.0
      Enddo

      If (tagAVERMOVE .EQ. 0) Then
      Do I = 1, NTOT3T - NPWALL
        X1(I) = X1(I) + U1(I) * DT
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = X2(I) + U2(I) * DT
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = X3(I) + U3(I) * DT
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
      Enddo
      Else ! tagAverMove
      Call AVER(U1, TEMP1)
      Call AVER(U2, TEMP2)
      Call AVER(U3, TEMP3)
      Do I = 1, NTOT3T - NPWALL
        X1(I) = X1(I) + TEMP1(I) * DT
        If (X1(I).GT.XMAX) X1(I) = X1(I) - CUBESIZE
        If (X1(I).LT.XMIN) X1(I) = X1(I) + CUBESIZE
        If (DIMM.EQ.DIMM3D) Then
          X2(I) = X2(I) + TEMP2(I) * DT
          If (X2(I).GT.XMAX) X2(I) = X2(I) - CUBESIZE
          If (X2(I).LT.XMIN) X2(I) = X2(I) + CUBESIZE
        Endif
        X3(I) = X3(I) + TEMP3(I) * DT
        If (X3(I).GT.XMAX) X3(I) = X3(I) - CUBESIZE
        If (X3(I).LT.XMIN) X3(I) = X3(I) + CUBESIZE
      Enddo
      Endif
      !Write(*,*) "Move", U1(1)

      If (tagRegroup .NE. 0) Then
        Call Regroup
      Endif

      Endif

      End Subroutine
